library(foreign)
library(arm)
library(tidyverse)

library(Hmisc)
library(lubridate)
library(broom)
library(transformr)
library(sandwich)
library(gtools)
library(jtools)
library(effects)
library(clubSandwich)
library(miceadds)
library(eurostat)
library(dplyr)
library(ggplot2)
library(stringr)
library(mapproj)
library(maps)
library(sf)
library(estimatr)
library(rnaturalearth)
library(rnaturalearthdata)

##### Maps
# getting geo-codes of all of Europe
world <- ne_countries(continent = 'europe', scale = "medium", returnclass = "sf")

# getting the region geo-codes from Eurostat
geo_spatial <- get_eurostat_geospatial(output_class = "sf", resolution = "60",
                                       nuts_level = "all", cache = TRUE, update_cache = FALSE,
                                       cache_dir = NULL)

# loading the MRP estimates
MRP <- read_csv("/Users/madselkjaer/Dropbox/WEALTHPOL_Research/Papers/WEP/replication files/MRP_estimates.csv")
MRP <- MRP %>%
  rename(NUTS_ID=X1)

# adding countries and regions that I want to be drawn but that are not in the ESS
m <- data.frame(matrix(NA, ncol = 3, nrow = 12))
names(m) <- c("NUTS_ID", "popu_index", "popu_index_standardized")
m$NUTS_ID <- c("LV", "RO", "HR", "AL", "SI", "RS", "ME", "MK", "ITF2", "LU", "LI", "TR")

# appending the data frames
MRP.m <- bind_rows(MRP,m)

# merging MRP estimates with geo codes
MRP_geo <- geo_spatial %>% 
  inner_join(MRP.m, by="NUTS_ID")

df.text.A <- data.frame(
  label = "A",
  x     = -16,
  y     = 70
)


# plotting the map
ggplot(data = world) +
  geom_sf(color = "dim grey", fill="white", size=.1)+
  geom_sf(data=MRP_geo, aes(fill=popu_index_standardized), color="dim grey", size=.1) +
  scale_fill_continuous(low="#f7fcf0", high = "#084081", name="", guide="colorbar", na.value="white") +
  labs(title="") +
  theme_light() + theme(legend.position="right") +
  coord_sf(xlim=c(-17,30), ylim=c(29,70)) +
  theme(title=element_text(size=8),
        axis.text=element_text(size=8),
        legend.position = "right",
        legend.text = element_text(size=8),
        legend.key.size = unit(2, 'lines'),
        axis.title.x=element_blank(), 
        axis.title.y=element_blank()) +
  guides(fill = guide_colorbar(barwidth = 0.25, barheight = 18))+
  geom_text(
    data    = df.text.A,
    mapping = aes(x = x, y = y, label = label), size=5, fontface="bold"
  )

### Workplace 
google <- readstata13::read.dta13("~/Dropbox/WEALTHPOL_Research/Papers/WEP/Data/Europe final data/ESS_GOOGLE.dta") 

google.dates <- google %>%
  filter(date_alt=="2020-02-20" | date_alt=="2020-04-30" | date_alt=="2020-06-25" )%>%
  rename(NUTS_ID=cregion)  

google.dates.geo <- geo_spatial %>% 
  inner_join(google.dates, by="NUTS_ID")


df.text.B <- data.frame(
  label = "B",
  date_alt   = "2020-02-20",
  x     = -16,
  y     = 70
)

# workplace
ggplot(data = world) +
  geom_sf(color = "dim grey", fill="white", size=.1)+
  geom_sf(data=google.dates.geo, aes(fill=workplace), color="dim grey", size=.1) +
  facet_wrap(.~date_alt, nrow=1)+
  scale_fill_continuous(low="#7f0000", high = "#fff7ec", guide="colorbar", name="") +
  labs(title="") +
  theme_light() +
  coord_sf(xlim=c(-17,30), ylim=c(29,70)) +
  theme(strip.background = element_rect(color="black", fill="white"),
        axis.text=element_text(size=8),
        axis.title=element_text(size=12),
        legend.position = "right",
        legend.text = element_text(size=8),
        strip.text.x = element_text(size = 8, color="black"),
        legend.key.size = unit(2, 'lines'),
        axis.title.x=element_blank(), 
        axis.title.y=element_blank())+
  guides(fill = guide_colorbar(barwidth = 0.25, barheight = 18))+  
  geom_text(
    data    = df.text.B,
    mapping = aes(x = x, y = y, label = label), size=5, fontface="bold"
  )



